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Given a quantum Hamiltonian and its evolution time, the corresponding unitary evolution op¬ 
erator can be constructed in many different ways, corresponding to different trajectories between 
the desired end-points and different series expansions. A choice among these possibilities can then 
be made to obtain the best computational complexity and control over errors. It is shown how a 
construction based on Grover’s algorithm scales linearly in time and logarithmically in the error 
bound, and is exponentially superior in error complexity to the scheme based on the straightfor¬ 
ward application of the Lie-Trotter formula. The strategy is then extended first to simulation of 
any Hamiltonian that is a linear combination of two projection operators, and then to any local 
efficiently computable Hamiltonian. The key feature is to construct an evolution in terms of the 
largest possible steps instead of taking small time steps. Reflection operations and Chebyshev ex¬ 
pansions are used to efficiently control the total error on the overall evolution, without worrying 
about discretisation errors for individual steps. We also use a digital implementation of quantum 
states that makes linear algebra operations rather simple to perform. 
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I. INTRODUCTION 

Richard Feynman advocated development of quantum 
computers as efficient simulators of physical quantum 
systems [l|. Real physical systems are often replaced by 
simplified models in order to understand their dynamics. 
Even then, exact solutions are frequently not available, 
and it has become commonplace to study the models us¬ 
ing elaborate computer simulations. Classical computer 
simulations of quantum models are not efficient—well- 
known examples range from the Hubbard model to lattice 
QCD—and Feynman argued that quantum simulations 
would do far better. 

As a concrete realisation of Feynman’s argument, it is 
convenient to look at Hamiltonian evolution of a many- 
body quantum system. Quantum simulations can sum 
multiple evolutionary paths contributing to a quantum 
process in superposition at one go, while classical sim¬ 
ulations need to evaluate these paths one by one. For¬ 
malisation of this advantage, in terms of computational 
complexity, has gradually improved over the years. Real 
physical systems are governed by local Hamiltonians, i.e. 
where each component interacts only with a limited num¬ 
ber of its neighbours independent of the overall size of 
the system. Lloyd constructed a quantum evolution al¬ 
gorithm for such systems Q, based on the discrete time 
Lie-Trotter decomposition of the unitary evolution oper¬ 
ator, and showed that it is efficient in the required time 
and space resources. Aharonov and Ta-Shma rephrased 
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the problem as quantum state generation, treating the 
terms in the Hamiltonian as black box oracles, and ex¬ 
tended the result to sparse Hamiltonians in graph the¬ 
oretical language Q- The time complexity of the al¬ 
gorithm was then improved jdj-jfij], using Suzuki’s higher 
order generalisations of the Lie-Trotter formula [7] and 
clever decompositions of the Hamiltonian. Recently the 
error complexity of the evolution has been reduced from 
power-law to logarithmic in the inverse error, using the 
strategy of discrete time simulation of multi-query prob¬ 
lems [8j. This is a significant jump in computational com¬ 
plexity improvement that needs elaboration and under¬ 
standing. In this article, we explicitly construct efficient 
evolution algorithms first for Hamiltonians that are lin¬ 
ear combinations of two projection operators, to expose 
the physical reasons behind the improvement , and then 
extend the strategy to any local efficiently computable 
Hamiltonian. Our constructive methods differ from the 
reductionist approach of Ref.jgj, we improve upon ear¬ 
lier results, and clearly demonstrate how the algorithms 
work in practice. 

Computational complexity of a problem is a measure 
of the resources needed to solve it. Conventionally, the 
computational complexity of a decision problem is spec¬ 
ified in terms of the size of its input, noting that the size 
of its output is only one bit. This framework is extended 
to problems with different output requirements (e.g. find 
the optimal route for the travelling salesman problem or 
evaluate n to a certain precision), by setting up succes¬ 
sive verifiable bounds on the outputs. For example, the 
problem of evaluating 7r can be implemented as first con¬ 
firming that 7r £ [3,4], and then narrowing down the in¬ 
terval by bisection, adding one bit of precision for every 
decision made. In such a scenario, the number of deci- 
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sion problems solved equals the number of output bits, 
and the complexity of the original problem is the sum of 
the complexities for the individual decision problems. It 
is therefore appropriate to specify the complexity of the 
original problem in terms of the size of its input as well 
as its output, especially for function evaluation problems. 
Generalizing the conventional classification, the compu¬ 
tational algorithm can then be labeled efficient if the re¬ 
quired resources are polynomial in terms of the size of 
both its input and its output. We label such algorithms 
as belonging to the class P:P, explicitly expressing their 
computational complexity with respect to their input as 
well as output sizes. Simultaneous consideration of both 
input and output dependence of complexity is natural 
for reversible computation. It is also necessary when ex¬ 
tending finite precision analog computation to arbitrary 
precision digital computation. Note that our definition 
of P:P efficiency differs from the concept of “simulatable 
Hamiltonians” in Ref. Q. 

The traditional computational complexity analysis 
(e.g. P vs. NP classification) does not discuss much how 
the complexity depends on the output precision, and the 
task is relegated to design of efficient methods for ar¬ 
bitrary precision numerical analysis [9]. We stress that 
both input and output size dependence of the compu¬ 
tational complexity are equally important for practical 
function evaluation problems. Popular importance sam¬ 
pling methods are not efficient according to our criterion, 
because the number of iterations needed in the compu¬ 
tational effort has a negative power-law dependence on 
the precision e (i.e. ./Viter oc e~ 2 as per the central limit 
theorem). On the other hand, finding zeroes of a func¬ 
tion by bisection is efficient (i.e. iViter oc loge), and 
finding them by Newton’s method is super-efficient (i.e. 
Wter OC log loge). 

We describe the Hamiltonian simulation problem in 
Section II, along with the important ingredients re¬ 
quired for its optimisation. In Section III, we formu¬ 
late the database search problem as Hamiltonian evolu¬ 
tion. While the evolution results are well-known (see for 
instance Ref. 0,0), we focus on the error complex¬ 
ity which has not been optimised in the literature. Our 
analysis explicitly shows how Grover’s large time-step al¬ 
gorithm is exponentially superior to small time-step al¬ 
gorithms approximating continuous time evolution. It 
also demonstrates that the large time-step algorithm ef¬ 
fectively simulates a very different Hamiltonian than the 
small time-step algorithm, but yields the same total evo¬ 
lution operator [12j. As an important ingredient, we in¬ 
troduce digital representation for quantum states that 
makes performing linear algebra operations with them 
straightforward. In Section IV, we construct a series ex¬ 
pansion evolution algorithm for Hamiltonians that are 
linear combinations of two projection operators. We 
carry out a partial summation of the series, and demon¬ 
strate how evaluation of a truncated series of large-step 
reflection operators improves the error complexity expo¬ 
nentially compared to the small-step Lie-Trotter formula. 


Our analytic results are supported by numerical tests. 
Finally in Section V, We combine the methods of Cheby- 
shev series expansion and digital representation, to con¬ 
struct an efficient simulation algorithm for any local ef¬ 
ficiently computable Hamiltonian. We conclude with an 
outlook for our methods, and some general results for 
projection operators are collected in an Appendix. 


II. QUANTUM HAMILTONIAN SIMULATION 


The Hamiltonian simulation problem is to evolve an 
initial quantum state |'i/>(0)) to a final quantum state 
\if{T)) 1 in presence of interactions specified by a Hamil¬ 
tonian H(t): 


m) 

U(T) 


u(T)\m >, 


exp 



H(t)dt ) 


(1) 


The initial state can often be prepared easily, while the fi¬ 
nal state is generally unknown. The path ordering of the 
unitary evolution operator U(T), denoted by the sym¬ 
bol V in Eq.dlJ, is necessary when various terms in the 
Hamiltonian do not commute. The properties of the final 
state are subsequently obtained from expectation values 
of various observables: 


(Oa) = {m\O a \m) . ( 2 ) 

In typical problems of quantum dynamics, both these 
parts—the final state and the expectation values—are 
determined probabilistically upto a specified tolerance 
level. They also require different techniques, and so it 
is convenient to deal with them separately. In this ar¬ 
ticle, we focus only on the former part; the latter part 
has been addressed in Refs. 0,0], and still needs expo¬ 
nential improvement in the dependence of computational 
complexity on the output precision to belong to the class 
P:P. For simplicity, we also restrict ourselves to problems 
where both the Hamiltonian H and the observables O a 
are bounded 0 ]. 

It is also possible to define the Hamiltonian simulation 
problem as the determination of the evolution operator 
U(T), and omit any mention of the initial and the final 
states. The accuracy of the simulation is then specified 
by the norm of the difference between simulated and ex¬ 
act evolution operators, say || U(T) — U(T)\\ < e. In 
actual implementation, the simulated U{T ) may not be 
exactly unitary, due to round-off and truncation errors, 
but the preceding measure for the accuracy of the simu¬ 
lation still suffices as long as e is small enough. 

We concern ourselves here only with Hamiltonians act¬ 
ing in finite /V-dimensional Hilbert spaces. A general 
Hamiltonian would then be a dense N x N matrix, and 
there is no efficient way to simulate it. So we restrict 
the Hamiltonian according to the following features com¬ 
monly present in physical problems: 
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(1) The Hilbert space is a tensor product of many small 
components, e.g. N = 2 n for a system of n qubits. 

(2) The components have only local interactions irrespec¬ 
tive of the size of the system, e.g. only nearest neighbour 
couplings. That makes the Hamiltonian sparse, with 
O(N) non-zero elements. 

(3) The Hamiltonian is specified in terms of a finite num¬ 
ber of efficiently computable functions, while the argu¬ 
ments of the functions can depend on the components, 
e.g. the interactions are translationally invariant. 

These features follow the notion of Kolmogorov complex¬ 
ity, where the computational resources needed to describe 
an object are quantified in terms of the compactness 
of the description. With a compact description of the 
Hamiltonian, the resources needed to just write it down 
do not influence the simulation complexity (Till . 

Such sparse Hamiltonians can be mapped to graphs 
with bounded degree d, with the vertices representing 
physical components of the system and the edges denot¬ 
ing the interactions between neighbouring components. 
Their simulations can be easily parallelised—on classical 
computers, Hamiltonians with these features allow SIMD 
simulations with domain decomposition. 

We note that long range physical interactions do ex¬ 
ist, but simulation of generic dense Hamiltonians is not 
efficient E3- Only with some extra properties, dense 
Hamiltonians can lead to non-local evolution operators 
having compact descriptions. A useful example is FFT, 
which describes a dense but factorisable unitary trans¬ 
formation that can be efficiently implemented, but we do 
not consider such possibilities here. 

With all these specifications, efficient Hamiltonian 
simulation algorithms in the class P:P use computational 
resources that are polynomial in log(lV), d and log(e). 


A. Hamiltonian Decomposition 

Efficient simulation strategy for Hamiltonian evolution 
has two major ingredients. In general, exponential of 
a sparse Hamiltonian is not sparse, which makes exact 
evaluation of exp (—iHt) difficult. So the first ingredient 
is to decompose the sparse Hamiltonian as a sum of non¬ 
commuting but block-diagonal Hermitian operators, i.e. 
H = Y^i =i The motivation for such a decomposition 

is twofold: 

(a) Functions of individual Hi, defined as power series, 
can be easily and exactly calculated for any time evolu¬ 
tion t, and they retain the same block-diagonal structure. 

(b) The blocks are decoupled and so can be evolved si¬ 
multaneously, in parallel (classically) or in superposition 
(quantum mechanically). 

Furthermore, the blocks can be reduced in size all the way 
to a mixture of 1 x 1 and 2x2 blocks. The lxl blocks 
just produce phases upon exponentiation, while the 2x2 
blocks can be expressed as linear combinations of identity 
and projection or reflection operators (i.e. (1 + h ■ <x)/2 or 
h ■ a respectively, where h is a unit vector and eq are the 


three Pauli matrices). There is no loss of generality in 
such a choice; it is just a convenient choice of basis that 
simplifies the subsequent algorithm. Projection or reflec¬ 
tion operators with only two distinct eigenvalues can be 
interpreted as binary query oracles. Their large spectral 
gaps also help in rapid convergence of series expansions 
involving them. 

In general, Hi can be systematically identified by an 
edge-colouring algorithm for graphs Q, with distinct 
colours (labeled by the index i) for overlapping edges. 
As per Vizing’s theorem, any simple graph of degree d 
can be efficiently coloured with d+1 colours. Physical 
models are often defined on bipartite graphs, for which 
the colouring algorithms are simpler than those for gen¬ 
eral graphs and need d colours. Identification of Hi also 
provides a compressed labeling scheme that can be used 
to address individual blocks. 

Actual calculations do not need explicit construction 
of U(T), rather only the effect of U(T) on the quantum 
state IV'(O)) has to be evaluated. That is accomplished by 
breaking down the calculation into steps, each of which 
consists of the product of a sparse matrix with a vec¬ 
tor, e.g. exp(— iHiffil'tp). The simulation complexity is 
then conveniently counted in units of such sparse matrix- 
vector products. 

As a simple illustration, the discretised Laplacian for 
a one-dimensional lattice has the block-diagonal decom¬ 
position given by: 
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( 3 ) 


This decomposition, H — H a + H e , has the projection 
operator structure following from H% = 2 H a and FT 2 = 
2 H e . Graphically, the break-up can be represented as: 


o e o e o 


where H a and H e are identified by the last bit of the 
position label. Eigenvalues of H are 4sin 2 (fc/2) in terms 
of the lattice momentum k , while those of H 0 and H e are 
just 0 and 2. 
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B. Evolution Optimisation 


Given that individual Hi can be exponentiated exactly 
and efficiently, their sum H can be approximately expo¬ 
nentiated using the discrete Lie-Trotter formula: 

exp ( — iHT) = exp HE HiT ) (4) 

i 

sa ^ exp (—iHiAt)^ , m = T/At. 


This replacement maintains unitarity of the evolution ex¬ 
actly, but may not preserve other properties such as the 
energy. The accuracy of the approximation is commonly 
improved by making At sufficiently small, sometimes ac¬ 
companied by higher order discretisations pjj]. This ap¬ 
proach has been used for classical parallel computer sim¬ 
ulations of quantum evolution problems [l9|, [20] . 

In contrast, the second ingredient of efficient Hamilto¬ 
nian simulation is to use as large At as possible. When 
the exponent is proportional to a projection operator, 
the largest At is the one that makes the exponential a 
reflection operator [2lJ . In general, use of any fixed con¬ 
stant At changes the leading scaling behaviour of the 
error complexity from a power-law dependence on e to a 
logarithmic one. The extreme strategy of choosing the 
largest possible At not only keeps the evolution accurate 
by reducing the round-off and the truncation errors, but 
also optimises the scaling proportionality constant |22| . 
It is not at all obvious how such a result may arise, and so 
we demonstrate it first in Section III using the database 
search problem as an explicit example, and then in Sec¬ 
tion IV for Hamiltonians that are a linear combination 
of two more general projection operators. 

Our efficient Hamiltonian simulation algorithms, ex¬ 
plicitly constructed using the two ingredients just de¬ 
scribed, have computational complexity 


O 


log jt/e) 

log(log(f/e)) 



(5) 


Here C is the computational cost of a single time step, 
which only weakly depends on t and e. It is C that char¬ 
acterises how computational complexity of classical im¬ 
plementation is improved in the quantum case, through 
conversion of independent parallel execution threads into 
quantum superposition. We point out that to keep the 
discretisation error under control, digital calculations 
need 6-bit precision, with b = Q(log((t/e) log(f/e))). For 
Lsparse Hamiltonians whose elements can be evaluated 
efficiently, the computational cost C is 0(lNb 3 ) classically 
and 0(lnb 3 ) quantum mechanically. 


identified with the individual items. It takes an initial 
state whose amplitudes are uniformly distributed over 
all the items, to the target state where all but one am¬ 
plitudes vanish. Let {|*}} be the set of basis vectors, |s) 
be the initial uniform superposition state, and 1 1) be the 
target state corresponding to the desired item. Then 

\m) = \s) , mn = 1 1) , 

|(i|s}| = 1/VN , (*|i) = S it . (6) 


The simplest evolution schemes taking |s) to 1 1) are 
governed by time-independent Hamiltonians that depend 
only on |s) and \t). The unitary evolution is then a rota¬ 
tion in the two-dimensional subspace, formed by |s) and 
|t), of the whole Hilbert space. In this subspace, let 


!<>=(;) • 

1 ' \y/{N-l)/Nj ■ 


(7) 


On the Bloch sphere representing the density matrix, the 
states |s)(s| and |f)(f| are respectively given by the unit 
vectors h s = ( 2v/ j^ -1 ,0, — l) and h t = (0,0,1). The 

angle between them is cos _1 (-^- — 1), which is twice the 
angle cos _1 (l/\/]V : ) between |s) and |t) in the Hilbert 
space. 

For a time-independent Hamiltonian, the time evolu¬ 
tion of the state is a rotation at a fixed rate around a 
direction specified by the Hamiltonian: 


U(t) = exp (—iHt) = exp(—in# • d c ot) . (8) 

For the database search problem, U(T)\s) = 11), upto 
a phase arising from a global additive constant in the 
Hamiltonian. There are many possible evolution routes 
from the initial to the target state, and we consider two 
particular cases in turn. 


A. Farhi-Gutmann’s and Grover’s Algorithms 

Grover based his algorithm on a physical intuition 
[23[, where the potential energy term in the Hamilto¬ 
nian attracts the wavefunction towards the target state 
and the kinetic energy term in the Hamiltonian diffuses 
the wavefunction over the whole Hilbert space. Both the 
potential energy |f)(f| and the kinetic energy |s)(s| [24| 
terms are projection operators. The corresponding time- 
independent Hamiltonian is 


III. QUANTUM DATABASE SEARCH AS 
HAMILTONIAN EVOLUTION 


H c = |s)(s| + \t){t\ 


( i + A AAEIN 

1 ^ N N 

I yiyrri , _ i_ I 
\ N N J 

T y/N~^l 1 

J + —Iv— 


(9) 


The quantum database search algorithm works in an 
V-dimensional Hilbert space, whose basis vectors are 


N 
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That gives rise to the evolution operator (omitting the 
global phase) 

Uc{t) = exp ( — in ■ a t/VN) 

= cos(t/y/N) — in ■ asin(t/y/N) , (10) 

which is a rotation by angle 2 t/y/N around the direction 
n=(y/(N-l)/N,0,l/VN) T on the Bloch sphere. 

The (unnormalised) eigenvectors of He are |s) ± |t), 
which are the orthogonal states left invariant by the evo¬ 
lution operator Uc(t). On the Bloch sphere, their den¬ 
sity matrices point in the directions ±h, which bisect the 
initial and the target states, |s)(s| and \t)(t\. Thus a ro¬ 
tation by angle it around the direction h takes |s)(s| to 
|i)(t| on the Bloch sphere. In the Hilbert space, the ro¬ 
tation angle taking |s) to 1 1) is then n/2, and so the time 
required for the Hamiltonian search is T — ( n/2)y/N [25| . 

Grover made an enlightened jump from this scenario, 
motivated by the Lie-Trotter formula. He exponentiated 
the projection operators in He to reflection operators; 
R = exp(±*7rP) = 1 — 2 P for any projection operator P. 
The optimal algorithm that Grover discovered iterates 
the discrete evolution operator [26} . 

U G = 


With Ug = exp(— iHar), it corresponds to the evolution 
Hamiltonian 


-(l-2|a)<a|)(l-2|t>(t|) 


1 - 


TV 


/jV-1 

TV 


o VN^l 2 

. Z TV 1 TV 


(ii) 


, 2 . „ y/N — 1 

(1 --)/ + 2 ? ——u. 


and the evolution step 


T 


N 

y/N^l 
2 N 

y/W^T 


sin 

sin 



(13) 


It is an important non-trivial fact that He is the com¬ 
mutator of the two projection operators in He'- 

ir G = *[|t)<t|,| s ><*|] . (14) 


This commutator is the leading correction to the Lie- 
Trotter formula in the Baker-Campbell-Hausdorff (BCH) 
expansion [27}, making Grover’s algorithm an ingenious 
summation of the BCH expansion for the evolution op¬ 
erator. 

On the Bloch sphere, each Ug step is a rotation by 
angle 2 ry/N — 1/N = 4sin^ 1 (l/v / TV) around the direc¬ 
tion fie = (0,1,0) T , taking the geodesic route from the 
initial to the final state. That makes the number of steps 
required for this discrete Hamiltonian search, 


Qt 


i cos ' 1 (|- 1 )/ sin "‘('^) 

cos-^l/y/N) _ tt ^ 

2sin- 1 (l /y/N) ~ 4 


(15) 


Note that h and uq are orthogonal, so the evolution 
trajectories produced by rotations around them are com¬ 
pletely different from each other, as illustrated in Fig|T| 
It is only after a specific evolution time, corresponding 
to the solution of the database search problem, that the 
two trajectories meet each other [28j- The evolution op¬ 
erators in the two cases are: 



FIG. 1: Evolution trajectories on the Bloch sphere for the 
quantum database search problem, going from |s) to |t). The 
Hamiltonians He and Hg generate rotations around the di¬ 
rections n and tig respectively. 


FIG. 2: Quantum logic circuit for the fractional query oracle 
operator O$ = exp(i</>|t)(t|). The oracle flips the ancilla bit 
iff its input is the target state, and the standard binary query 
oracle operator corresponds to <j> = it. 
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The rates of different Hamiltonian evolutions can be 
compared only after finding a common convention to fix 
the magnitude and the ease of simulation of the Hamil¬ 
tonians: 

(1) Magnitude: A global additive shift of the Hamilto¬ 

nian has no practical consequence. So possible compari¬ 
son criteria can be the norm of the traceless part of the 
Hamiltonian or the spectral gap over the ground state. 
H-ffcll = 1 /y/N and ||.Hg|| = are comparable, with 

the same limit as N —> oo. 

(2) Ease of simulation: The binary query oracle can be 
used to produce various functions of \t){t\. 

(a) He can be easily simulated by alternating small 
evolution steps governed by |s)(sjand |t)(t|, accord¬ 
ing to the Lie-Trotter formula [29|. Each evolution 
step governed by \t)(t\ needs two binary query ora¬ 
cles, as shown in FigJ^l 

(b) Ug is easily obtained using one binary query 
oracle per evolution step. 


B. Equivalent Hamiltonian Evolutions 

When different evolution Hamiltonians exist, corre¬ 
sponding to different evolution routes from the initial to 
the final states, one can select an optimal one from them 
based on their computational complexity and stability 
property. This feature can be used to simplify the Hamil¬ 
tonian evolution problem by replacing the given Hamil¬ 
tonian by a simpler equivalent one. Two Hamiltonian 
evolutions are truly equivalent, when their correspond¬ 
ing unitary evolution operators are the same (for a fixed 
evolution time and upto a global phase). The intersec¬ 
tion of the two evolution trajectories is then independent 
of the specific initial and final states. 

For the database search problem, we observe that 

U c (T) = i(l-2\t)(t\) (U g ) Qt ■ (18) 

So with an additional binary query oracle, one evolu¬ 
tion can be used as an alternative for the other, with¬ 
out worrying about the specific choices of |s) and 1 1). 
(The additional oracle is needed to make the two Hamil¬ 
tonian evolutions match, although it is not required for 
the database search problem.) 

For a more general evolution time 0 < t < T, we have 
the relation (analogous to Euler angle decomposition), 

Uc(t) = exp (i/3cr 3 ) (Uc) Qt ex P (*(f + , (19) 


Since 03 = 2|i)(i| — 1, each phase rotation is a fractional 
query oracle and can be obtained using two oracle calls 
[291 ] . The parameters in Eq. (flUl) are given by 


Qt — 


sin t 


2sin” 1 (l /y/N) 


/3 = — — — - tan 


2 ’ 
t 


(7v lan 7v)' (20) 


They yield 



In — i . t \ \ 

( i sm 1 

V N Sm VNr 2 ) 


/C0S2 VN + N Sm2 Jn 


( 21 ) 


N n Xsm Jn 


s i n _L_ 

n sm JN 


cos2 Jn + n sm “ Jn . 


whose elements are the same as those of Uc{t) upto phase 
factors. 

Thus Hq can be used to obtain the same evolution as 
He, even though the two Hamiltonians are entirely dif¬ 
ferent in terms of their eigenvectors and eigenvalues — 
a rare physical coincidence indeed! A straightforward 
conversion scheme is to break up the duration of evolu¬ 
tion for He into units that individually solve a database 
search problem, simulate each integral unit according to 
Ea. (fl 8 l) . and the remaining fractional part according to 
Eq.®. 


C. Unequal Magnitude Evolution Operators 


Now consider the generalisation of He to the situation 
where the coefficients of |s)(s| and \t)(t\ are unequal. In 
that case, the rotation axis for continuous time evolution 
is not the bisector of the initial and the target states, even 
though it remains in the oq — 03 plane. As a result, one 
cannot reach the target state exactly at any time. The 
database search succeeds only with probability less than 
one, although the rotation angle on the Bloch sphere for 
the closest approach to the target state remains n. The 
equal coefficient case is therefore the choice that max¬ 
imises the database search success probability. 

On the other hand, the results obtained using discrete 
time evolution for the Hamiltonian simulation problem 
are easily extended to situation where the two projec¬ 
tion operators have unequal coefficients. Without loss of 
generality, we can choose 


H 


a|s)(s| + |t)<t| 



ay/N^l 

N 


o\ + 



( 22 ) 

fh- 


i.e. Uc(t) can be generated as Q t iterations of the Grover with real a £ [—1,1]. It gives rise to the evolution oper- 

operator Ug, preceded and followed by phase rotations. ator (without the global phase) 
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U(t) 



i (^ + 77) sin {At) 
EvE sm(At) 


where 


A 2 



(24) 



N — 2 \2 

~^r a ) + 


(N - 1 )a 2 
N 2 


> 


sa\/N — 1\ 2 

V N ) 


The rotation axis for this evolution is still in the £ 71-03 
plane, and the commutator of the two terms in the Hamil¬ 
tonian is still proportional to Hq- As a consequence, U (t) 
can still be expressed as Q iterations of the Grover oper¬ 
ator Ug, preceded and followed by phase rotations: 

U(t) = exp (i/3a 3 ) ( U G ) Q exp + /3)<t 3 ^J , (25) 


with the parameters given by 

Q = 8111-1 ***>)/ (2sln ~ 1 ’ 

P = tan ~ l ((“ 2 “ + ^ tan(At)) .(26) 

Note that for Q < 0, we need to iterate the operator 
U^ = ^l^2\t)mi~2\s){s\). 


D. Discretised Hamiltonian Evolution Complexity 

In a digital implementation, all continuous variables 
are discretised. That allows fault-tolerant computation 
with control over bounded errors. But it also introduces 
discretisation errors that must be kept within specified 
tolerance level by suitable choices of discretisation inter¬ 
vals. When Hamiltonian evolution is discretised in time 
using the Lie-Trotter formula, the algorithmic error de¬ 
pends on At, which has to be chosen so as to satisfy the 
total error bound e on U(t). The overall computational 
complexity is then expressed as a function of t and e. 

For the simplest discretisation, 

1 

exp ( - * E E ' At ) 

1=1 

= exp ( — iHiAt) ... exp ( — iHiAt) (27) 
x exp ( — iE^ 2 \At) 2 ), 

E {2) = i -J2 ^ H >\+°( At )- ( 28 ) 

i<j 


cos(Ai) + i(ij! + *) s in(ylt)y 


For the symmetric discretisation, 

1 

exp ^ — i E HiAt ) 

2=1 

= ^ exp (—iHi At/2)... exp(— iH 1 At/2)^ 
x ^ exp(— iHiAt/2)... exp (—iHi At/2)^ (29) 

x exp ( — iE^ (At) 3 ) , 

Ei3) = + 
i<j 

+ ^ E (2 [H i ,[H J ,H k ]] + [H j ,[H i ,H k ]]) 

i<j<k 

+ O(At) . (30) 

Here quantify the size of the discretisation error. 
These discretisations maintain exact unitarity, but do not 
preserve the energy when H and E^ do not commute. 

For any unitary operator X, the norm ||X|| is equal to 
one (measured using either Tr(X^X) or the magnitude 
of the largest eigenvalue). That makes, using Cauchy- 
Schwarz and triangle inequalities, 

l|x m -y m || = ||(x-y)(x m - 1 + ... + y m " 1 )|| 

< m\\X-Y\\ . (31) 

So for the total evolution to remain within the error 
bound ei, we need 

m\\exp(-iE {k \At) k )- I\\ ~ m \\E w \\(At) k (32) 

= m l - k t k \\E (k) \\ < ci . 

The error probability can be rapidly reduced by re¬ 
peating the evolution a multiple number of times, and 
then selecting the final result by the majority rule (not 
as average). This simple procedure produces an error 
bound similar to higher order discretisation formulae. 
With R repetitions, the error probability becomes less 
than which can be made smaller than any 

prescribed error bound e [30j |. With exact exponentia¬ 
tion of the individual terms Hi , the computational cost 
to evolve for a single time step At, i.e. C, does not de¬ 
pend on At. Thus the complexity of the Hamiltonian 
evolution becomes 

o(mRC) = l ^X!mn RC ) ■ < 33 > 
















With superlinear scaling in t and power-law scaling in 
e, this scheme based on small At is not efficient. Note 
that for the Hamiltonian He, l = 2 is fixed, and both 
||£( 2 )|| and \\E^\\ are C^TV" 1 / 2 ). So for evolution 
time T = ©(TV 1 / 2 ), the time complexity becomes lin¬ 
ear, 0{Te~ 1 /^ k ~ 1 ' , ^ R / 2 ^RC), while power-law scaling in 
e remains unchanged. 

Grover’s optimal algorithm uses a discretisation for¬ 
mula where exp(—TiL^Af) are reflection operators. The 
corresponding time step is large, i.e. A tG = n for Eq.((4| 
applied to Eq.©. The large time step introduces an¬ 
other error because one may jump across the target state 
during evolution instead of reaching it exactly. Q t is not 
an integer as defined in Eq. (12(71) , and needs to be replaced 
by its nearest integer approximation [Qt + ^J in practice. 
For instance, the number of time steps needed to reach 
the target state in the database search problem is 



Since each time step provides a rotation by angle a = 
2sin~ 1 (l/\/TV) along the geodesic in the Hilbert space, 
and one may miss the target state by at most half a ro¬ 
tation step, the error probability of Grover’s algorithm is 
bounded by sin 2 (a/2) = 1/N. Since the preceding and 
following phase rotations in Ea. (fT9]l are unitary opera¬ 
tions, this error bound applies to Uc(t) as well. Once 
again, reducing the error probability with R repetitions 
of the evolution and the majority rule selection, we need 
2 r ~ 1 /N^ r / 2 ^ < e. The computational complexity of the 
evolution is thus 


o(Q t RCo) = o(i(-Af) Co ) 

= (35) 

With linear scaling in time and logarithmic scaling in e, 
this algorithm is efficient. 

It is easy to see why the two algorithms scale rather 
differently as a function of e. The straightforward appli¬ 
cation of the Lie-Trotter formula makes the time step At 
depend on e as a power-law. The total error of the algo¬ 
rithm is proportional to the total number of time steps 
to, and the resultant computational complexity then has 
a power-law dependence on e. The power can be reduced 
by higher order discretisations or by multiple evolution¬ 
ary runs and the majority rule selection, but it cannot 
be eliminated. On the other hand, with a large time step 
that does not depend on e, Grover’s algorithm has an er¬ 
ror that is independent of the evolution time. This error 
is easily suppressed by multiple evolutionary runs and 
the majority rule selection. The overall computational 
complexity is proportional to the number of evolutionary 
runs, which depends only logarithmically on e. 


E. Digital State Implementation 

To estimate the computational cost C, we need to spec¬ 
ify quantum implementation of linear algebra operations 
involving the block-diagonal operators Hi. It is routine to 
represent a quantum state in an TV-dimensional Hilbert 
space as 


jV-1 


N -1 


c ) = x i\fi ’ J2 


= i 


(36) 


3=0 


3=0 


where Xj are continuous complex variables. This analog 
representation is not convenient for high precision calcu¬ 
lations, and so we use the digital representation instead 


3l|, specified by the map 


3=0 


(37) 


This is a quantum state in a (2 h TV)-dimensional Hilbert 
space, where | Xj)b are the basis vectors of a 6-bit reg¬ 
ister representing the truncated value of Xj (a complex 
number Xj can be represented by a pair of real num¬ 
bers, and 2 b Xj are truncated to integers). This repre¬ 
sentation is fully entangled between the component in¬ 
dex state | j) and the register value state | Xj)b, with a 
unique non-vanishing \xj)t, (out of 2 b possibilities) for 
every \j). It is important to observe that no constraint is 
necessary on the register values in this representation— 
the perfect entanglement ensures unitary evolution in the 
(2 b TV)-dimensional space. This freedom allows simple im¬ 
plementation of linear algebra operations on \xj)i ,, trans¬ 
forming them among the 2 b basis states using only C-not 
and Toffoli gates of classical reversible logic, with the in¬ 
dex state | j) acting as control. For example, 


c\x) 


1 

Vn 


N-l 

b')l CX o)b 

3=0 


I x) + \y) —t 


1 

Vn 


N-l 

V2 \j)\ X 3+V3)b , 
3=0 


(38) 


map non-unitary operations on the left to unitary opera¬ 
tions on the right. The circuits described later in Section 
IIVB II and Section IV Cl combine such elementary opera¬ 
tions to construct power series. Note that a crucial re¬ 
quirement for implementing linear algebra operations in 
the digital representation is that only a single index ( “j” 
in the preceding formulae) controls the whole entangled 
state. 

The freedom to choose a convenient representation for 
the quantum states is particularly useful due to the fact 
that the quantum states are never physically observed. 
All physically observed quantities are the expectation 
values of the form in Eq.([2). So to complete the dig¬ 
ital representation, we need to construct for every ob¬ 
servable O a in the TV-dimensional Hilbert space a related 
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observable O a in the (2 b ./V)-dimensional Hilbert space, 
such that 

JV-l 

(x\O a \x) = ^ x j x l(j\°a\l) 

3,1 =o 

1 N ~ x 

= JjYl b(xj\(j\O a \l)\xi) b . (39) 

3,1=0 

For this equality to hold, it suffices to construct the op¬ 
erator O a = O 0 ® Ob, where the Hermitian operator Ob 
in the 2 b -dimensional Hilbert space satisfies 

(xj\Ob\xi) = Nx*xi . (40) 

For a single bit Xj , the solution is easily found to be the 
measurement operator O b= i = iV( 1 ~°’ 3 ). More gener¬ 
ally, we note that 

(zjKl + cri)® b |*i) = 1 , (41) 

and the place-value operator for a bit string, 

b-i _ 

V = Y^ 2 ~ k I m ® ( °~ 3 ) ® , (42) 

k=0 

gives V\xj) = Xj\xj). The solution to Ea. (l40l) . therefore, 
has a bit-wise fully factorised fornn independent of the 
quantum state and the observable [321 ]. 

O b = NV\l + a 1 )® b V . (43) 

The computational complexity of measurement of physi¬ 
cal observables in the digital representation is thus 0(b 2 ) 
times that in the analog representation. The advantages 
of arbitrary precision calculations and simple linear al¬ 
gebra, however, unambiguously favour the digital repre¬ 
sentation over the analog one. 

To efficiently incorporate the digital representation in 
the Hamiltonian simulation algorithm, methods must be 
found to not only manipulate the register values | Xj) effi¬ 
ciently, but also to initialise and to observe them. At the 
start of the calculation, we need to assume that the initial 
values Xj(0) can be efficiently computed from j. Then the 
initial state can be created easily using Hadamard and 
control operations, for N = 2", as 

1 N—l 

|0)|0)b H -= ]T |j)|0) b (44) 

v j=o 

N -1 

-j='E\j)\x j {0)h- (45) 

V 3=0 

When N is not a power of 2, some extra work is needed. 
A simple fix is to enlarge the j-register to the closest 
power of 2 and initialise the additional Xj to zero. There¬ 
after, the linear algebra operations can be implemented 
such that the additional Xj remain zero, and the overall 


normalisation (i.e. 1 /y/N) can be corrected in the final 
result as a proportionality constant. At the end of the 
calculation, we need to assume that the final state ob¬ 
servables, Eq.@, are efficiently computable from Xj(T). 
In such a case, the advantage of the digital representation 
is, as pointed out earlier, that the index j can be han¬ 
dled in parallel (classically) or in superposition (quantum 
mechanically). 

Digital computation with finite register size produces 
round-off errors, because real values are replaced by inte¬ 
ger approximations. To complete the analysis, we point 
out the standard cost estimate to control these errors. 
With 6 -bit registers, the available precision is 6 = 2 -b . 
Using simple-minded counting, elementary bit-level com¬ 
putational resources required for additions, multiplica¬ 
tions and polynomial evaluations are 0 ( 6 ), 0 ( 6 2 ) and 
0(6 3 ) respectively. (Overflow/underflow limit the de¬ 
gree of the polynomial to be at most 6 .) All efficiently 
computable functions can be approximated by accurate 
polynomials, so the effort needed to evaluate individual 
elements of Hi is 0 ( 6 3 ). 

The block-diagonal Hi can be exponentiated exactly. 
(Depending on the available quantum logic hardware, 
Euler angle decomposition may be used to convert rota¬ 
tions about arbitrary axes to rotations about fixed axes.) 
With fixed block sizes, exponential of any block of any Hi 
can therefore be obtained to 6 -bit precision with 0 ( 6 3 ) 
effort. The number of blocks is O(N), and so the clas¬ 
sical cost of multiplying exponential of Hi with a state 
is proportional to N. With an efficient labeling scheme 
for the blocks, the index j can be broken down into 0(n) 
tensor product factors (analogous to Ea. (l44l) l. and then 
quantum superposition makes the cost of multiplying ex¬ 
ponential of Hi with a state proportional to n. Thus 
the computational cost of an efficiently encoded matrix- 
vector product reduces from its classical scaling 0(Nb 3 ) 
to its quantum scaling 0(nb 3 ). 

For the database search problem, the number of ex¬ 
ponentiations of Hi needed for the Lie-Trotter formula is 
m(k—l)l, which reduces to 2 Q t for the Grover version. So 
with the choice m{k — 1)16 = O(e), i.e. 6 = fi(log(m/e)), 
the round-off error can be always made negligible com¬ 
pared to the discretisation error. The computational cost 
of a single evolution step then scales as 

C = 0(log V (log(f/e)) 3 ) = Cg , (46) 

and the overall evolution complexity of Eq. (1351) becomes 
0(-t loge (log(t/e)) 3 ). 


IV. FROM DATABASE SEARCH TO MORE 
GENERAL PROJECTIONS 

The Hamiltonian He for the database search prob¬ 
lem is a sum of two one-dimensional projection operators 
with equal magnitude. We next construct an accurate 
large time step evolution algorithm for the case where the 
two projection operators making up the Hamiltonian are 
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more than one-dimensional but block-diagonal, for ex¬ 
ample as in Eq.®. Our strategy now relies on a rapidly 
con verging series expansion, similar to the proposal of 
Ref. [33|, [3J] , instead of an equivalent Hamiltonian evolu¬ 
tion. We consider, in turn, series expansions in terms of 
projection operators and in terms of reflection operators. 


A. Algorithm with Projection Operators 

1. Series Expansion 

Consider the Hamiltonian decomposition 


at sufficiently high order. Also, the series in Eo. (l48l) can 
be efficiently summed using nested products, e.g. 




k =1 


= ffi(df + H 2 {c 2 I + ffi(c 3 J + ■ • ■))) ■ (52) 


The series in Ea. (148l) can be converted to a form related 
to the BCH expansion as: 


e ~iHt _ e —iH 1 t e —iH 2 t 


k =2 


4 2) (t)(fr 2 H 1 fl- 2 ...)fc) 


, (53) 


H = Hi + H 2 , Hi = H x , H\ = H 2 . (47) 

Then standard Taylor series expansion around t = 0 
yields 

OO 

exp (-iHt) = I+ ^c k (t) [ (H 1 H 2 H 1 ...) k 
fe =i 

+ {H 2 H x H 2 .. .)*] ,(48) 

with Chit = 0) = 0. Here {HiH 2 H\ ...)/. denotes a prod¬ 
uct of k alternating factors of Hi, starting with Hi. All 
other products of Hi reduce to the two terms retained on 
r.h.s. in Eo. (H51) . due to the projection operator nature 
of Hi, and the two have the same coefficients Cfc(f) by 
symmetry. It is worthwhile to observe that the structure 
of Ea. (1481) effectively sums up infinite series of terms— 
when truncated to order p , the series has 2p + 1 terms, 
compared to 2 P+1 — 1 terms in the corresponding series 
of Ref. 0. 

Differentiating Eq. (H51) . we obtain 
-i{Hi + H 2 ) x 

OO 

I + J2 Ck (*) ...) k + (H 2 H 1 H 2 .. .)*] 

k=l 

= E ^7^ {{HiHzHi . . ,) fc + {HzHiH 2 . . .) fc ] ,(49) 

k= 1 

which provides the recurrence relation for the coefficients, 

dc ^ = -i(c k (t ) + Cfc_i(t)) , (50) 

at 

with the initial condition Cq = 1. Iterative solution gives 


Noting that exp (iHit) = I + ( e lt — 1 )Hi, we can evaluate 
the coefficients as: 

c fc 1} = e * 4 ( c fc-i + Cfc) - c k -1 , 

cf = (e l ‘ - !) 2 (cfe -2 + c k - 1 ) + 4 1} . (54) 

Although the series in Ea. (l53l) starts with k = 2, c^\t) 
do not converge any faster than Ck{t) for larger k, and so 
there is no particular advantage in using it compared to 
the series in Eq. diSl) (35| . 

2. Complexity Analysis and Series Order Determination 

The computational complexity of Hamiltonian evo¬ 
lution using the Lie-Trotter formula is 0(mC), as in 
Ea. (133D . where m = t/At and C represents the computa¬ 
tional cost to evolve for a single time step. In order to 
keep the total evolution within the error bound e, At has 
to scale as a power of e, which in turn makes the com¬ 
putational complexity inefficiently scale as a power of e. 
Instead, with a truncated series expansion of Eo. (l48l) . we 
can choose At = 0(1). The order of series truncation, p, 
is then determined by the error bound e. A single time 
step needs 2 p nested linear algebra operations, with each 
operation consisting of a sparse matrix-vector product in¬ 
volving Hi, multiplication of a vector by a constant and 
addition of two vectors. So the computational complex¬ 
ity is 0(2mpC) with C denoting the computational cost 
of evaluating Hi and performing the linear algebra oper¬ 
ation. A simple quantum logic circuit to implement the 
linear algebra operation, using the digital representation 
of Section IIIIE1 is described later in Section IIVB II 
The series truncation error for a single time step, with 
1117,11 < 1 for projection operators, is 


Ck(t) = {-i) k e lt dtk ... / dt 2 / dti e l 

Jo Jo Jo 


= (- 1 ) 


k „—it 


2^ n I 

0=0 J ■ 


(51) 


Clearly \ck(t)\ = 0(t k /k\), and we can get as accurate 
approximations to e~ lHt as desired by truncating Eq. lBSl) 


A[exp(—iff At)] < 2 ^ |cfc(At)| . (55) 

/c—p+1 

It has to be bounded by e/m according to the triangle 
inequality. From Eq. (151l) . we have 


|c fc (Af)| = 


E 

j—k 


(iAty 


j'- 


< 


(A tf 
k\ 


1 - 


At 
k + 1 


■ (56) 
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The constraint deciding the order of series truncation is, 
therefore, 


2 TO 


(A t)P +1 
(P+ 1)! 



-2 

< e . 


(57) 


With At = 0(1), we have m = 0(t), and the formal solu¬ 
tion isp= 0(log(t/e)/log(log(t/e))). The computational 
complexity of the evolution is then 


0(2mpC) = O ( t 


log(*A) 

log(log(t/e)) 


(58) 


which makes the series expansion algorithm efficient. 

Finally, with block-diagonal Hi and finite precision cal¬ 
culations using 6-bit registers, the computational cost 
C is 0(nb 3 ). Then the choice mpS = O(e), i.e. b = 
12(log (mp/e)) = 12(log((t/e) log(t/e))), makes the round¬ 
off errors negligible compared to the truncation error. 


3. Unequal Magnitude Operators 

The series expansion algorithm is easily extended to 
the situation where the two projection operators appear¬ 
ing in the Hamiltonian have unequal coefficients. With 
H = a\H\ + a 2 H 2 , the series expansion takes the form 


B. Algorithm with Reflection Operators 

1. Series Expansion 


The series expansion can also be carried out in terms 
of the reflection operators Ri = I — 2 Hi , instead of the 
projection operators Hi. We then have 

e lt exp(-iHt) = exp (i(Ri + i? 2 )^) (63) 

OO 

= r 0 (t) I + ^2 r k(t) [(R1R2R1 •.■)* + (R2R1R2 • • ■)*] , 

k= 1 

with ro(t = 0) = 1 ,r k (t = 0) = 0. The structure of the 
terms in this series, with alternating reflection operators, 
is reminiscent of Grover’s algorithm. Differentiating this 
expansion, we obtain 


1 -{R 1 +R 2 )x (64) 

OO 

r 0 (t) 1 + ^2 r *;( 2 ) [(R1R2R1 ■■■)k + (R2R1R2 ■ ■ -)fc] 

fc=1 


dr 0 (t) ■A dr k (t) 

dt ^ dt 

fc =1 


[(R1R2R1 ■ ■ -)k + {R2R1R2 ■ ■ -)fe] 


exp (-iHt) = I + ^2 [ c k (t)(H 1 H 2 H 1 .. .) k 
k =1 

+ d k {t){H 2 HiH 2 ...) k ] ,(59) 

where c k {t = 0) = 0 = d k {t = 0), and c k (t) = d k {t) for 
even k by symmetry. Without loss of generality, one may 
choose ai € [—1,1], a 2 = 1 as in Ea. (l22l) . 

Differentiation of Eg. (l5f)l) leads to the recurrence rela¬ 
tions, 


which provides the recurrence relations for the coeffi¬ 
cients, 

^ = iriW ’^r = K rfc+lW+rfc - lW ) ’ (65) 

for k > 1. These are the recurrence relations for the 
Bessel functions. 

Explicit evaluation for the coefficient of identity in the 
series gives, using Rf = J, 


dc k (t ) 
dt 

dd k {t) 

dt 


= —ia\(c k (t) + dfe_i(t)) , 

= -ia 2 {c k -i(t) + d k (t)) , (60) 


with the initial conditions Co = 1 = do- These can be 
integrated to 

Ck (t) = -ia 1 e~' iait [ e iait 'd k -i(t')dt' , 
dk(t) = —ia 2 e~ ia2t f e ia2t c k -i(t')dt' , (61) 

and iteratively evaluated to any desired order. In partic¬ 
ular, for even fc, 


c k (t) = d k (t) = (a 2 c fc _i - ai<4_i)/(ai - a 2 ) . (62) 

With rapidly decreasing coefficients, |cfc(2)| = 0(t k /k\) = 
\d k (t)\, accurate and efficient truncations of Eg. (1551) are 
easily obtained. 



Thereafter, the recurrence relations determine r k [t) = 
i k J k {t). With \r k (t)\ = 0(t k /(2 k k\)), Eo. (l63l) converges 
significantly faster than Ea. (l48l) . It can also be summed 
efficiently using nested products. Furthermore, reflec¬ 
tions are unitary operators, and so they are easier to im¬ 
plement in quantum circuits than projection operators. 
These properties make Eq. flfldl) better to use in practice 
than Eo. flSl) . 

Summation of the series in Eo. (l63l) . truncated to order 
p , by nested products requires 2 p executions of the ele¬ 
mentary linear algebra operation fragment (r/ + i?)|a;). 
Each fragment contains three simple components: multi¬ 
plication of a vector by a unitary matrix, multiplication 
of a vector by a constant, and addition of two vectors. 
Its evaluation using the digital representation of Section 
1III El is schematically illustrated in FigJ31 Multiplication 
of a vector by a diagonal matrix, and addition of two 
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vectors are easy tasks. Multiplication of |a;) by the off- 
diagonal elements of Ri needs a little care, and can be ac¬ 
complished by shuffling the elements of |x). Since Ri are 
block-diagonal, this shuffling is only within each block, 
and requires a fixed number of permutations that depend 
on the block size but not on the system size. The time 
complexity of the series summation is thus O(pC), where 
C = 0(nb 3 ) using quantum superposition over the index 
j. The space resources required to combine together the 
results of all the fragments are a fixed number of n-bit 
registers and O(p) 6-bit registers. (The registers used 
for off-diagonal matrix multiplication in individual frag¬ 
ments can be reversibly restored to zero, and then reused 
in subsequent steps.) These features make the algorithm 
efficient, and the procedure is considerably simpler than 
the corresponding series summation method in Ref. (34] . 


< (At) p+1 ( _ At \ 

- 2 p+1 (p + 1)! \ 2(p + 2)J 


With t — rriAt, the order of series truncation is therefore 
decided by the constraint 


„ (Af) p+1 / At 

2 p +i(p+l)! V ~ 2(p + 2) 


-l 

< e . 


(70) 


For time step At = 0(1), the formal solution is again 
p = 0(log(i/e)/log(log(i/e))). That keeps the algorithm 
efficient, with the same computational complexity as in 

Eq.®. 


3. Unequal Magnitude Operators 


2. Complexity Analysis and Series Order Determination 

When the series of Eq. (l63l) is truncated at order p, with 
time step At and 11 Ri \ | = 1, the truncation error is 

OO 

A[exp(— iH At)] < 2 0^ |r*,(At)| . (67) 

k—p-\-1 

Since the Bessel functions obey 


s=0 v x 


it follows that (assuming (At) 2 < 8(p + 2)) 

(A t) k 


E i J *( A *)i ^ E 


k—p -\-1 


k—p -\-1 


2 k k\ 


( 68 ) 

(69) 


|r) 6 

10 

\j) 

10) 

|0)i 

|0)i 

I 0) b 

| 0>6 


FIG. 3: Digital quantum logic circuit for the linear algebra 
fragment | y) = (ri + Ri)\x) occurring in the nested evaluation 

of the series in Eg. (1631) . Among the controlled logic gates, H, 

I R I and 0 denote oracle operations specified by the Hamilto¬ 
nian and the initial state, while [X] stands for the generalised 
Toffoli gate implementing \a,b, c) —> \a,b, c + ab). The circuit 
is used with a uniform superposition over the index j. 
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When the two reflection operators have unequal coef¬ 
ficients in the Hamiltonian, we can expand 

exp (i(aiRi + a 2 i? 2 )^) = Po(t) I (71) 

OO 

+ E \Pk(t)(RlR2Rl •■•)*; + Qk(t){R2RlR2 ■ ■ -)fc] J 

fc=1 

with p 0 (t = 0) = 1 ,Pk(t = 0) = 0 = qk(t = 0), and 
Pfe(t) = qk(t) for even k by symmetry. 

Differentiation of Eq. m leads to the recurrence rela¬ 
tions, 


dpk(t) if ,, 

= - [a\qk-i(t) + a 2 qk+v 

= l -{a 2 Pk-i(t) + axpk+i 


dt 

dqk(t) 

dt 


(72) 


for k > 1. These can be iteratively solved to obtain the 
coefficients pk and qk for k > 2, to any desired accuracy, 
starting from the initial coefficients po = qo, p\ and q±. 
Explicit evaluation of these initial coefficients gives 


j =o y \z=o 


it \ 1^ 2 


2 ( 3-0 „ 2 1 


(73) 


Pi 


(*) = E 


i 


it \ 2 i+i 


3=0 


(2 j + 1) 




E 

o 


A (j + A 2(3-0+! 


I 


Jll 


(74) 


qi{t) is obtained from p\{t) by interchanging ai +>■ a 2 , 
and we also have the relation: 


dpoit) 

dt 


= ^oipi(t) + a 2 qi 


(75) 


The bounds |pfc(t)| = 0(t k /{2 k k\)) = |g fc (t)| make accu¬ 
rate and efficient truncations of Eq. (ED possible. 
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C. Numerical Tests 

The computational complexity bounds, Eg. (1551) and 
Eg. (1551) . have been obtained assuming that the evolu¬ 
tion errors during different time steps are unrelated. In 
practice, these bounds are not tight because correlations 
exist between evolution errors at successive time steps. 
To judge the tightness of the bounds, and also to esti¬ 
mate the scaling coefficients involved, we simulated the 
Lie-Trotter (with k = 2) and series expansion algorithms, 
Eos. (l27l) and (|48I63D respectively, for the one-dimensional 
discretised Laplacian ( H e + H a )/2 defined as per Eq.([5]). 

We carried out our tests on a one-dimensional periodic 
lattice of length L = 128, with a random initial state 
|^(0)}. We quantified the error as the norm of the dif¬ 
ference between the simulated and the exact states, i.e. 
e = |||^) — |V , )||- We also needed mp5 < e to keep the 
round-off errors under control. That was not possible 
with 32-bit arithmetic, and we used 64-bit arithmetic. 

We selected t = 100 to study the dependence of the 
error on the evolution step size and the series truncation 
order. Our results are displayed in Fig{4j For the series 
expansion algorithms, as expected, we observe that (a) 
the truncation order p depends linearly on log(e), (b) 
the numerical values are consistent with the bounds in 
Eos. (157l70D but the bounds are not very tight, and (c) 
the reflection operator series converges faster than the 
projection operator series. For At = 1, the coefficients 
Cfe(At) and rk(A t) decrease monotonically, and the series 
reach a given error e with a smaller order p compared 
to the case At = i r. But in the overall computational 



FIG. 4: Dependence of the error e on the truncation order p 
for the series expansion algorithms, and log 2 m for the Lie- 
Trotter algorithm. The symbols □, + and 0 respectively 
represent the results for the reflection operator series, the pro¬ 
jection operator series and the Lie-Trotter algorithm. Con¬ 
tinuous and dashed lines connect series expansion results for 
At = 1 and At = n respectively. 


complexity, this reduction in p (roughly a factor of 1.6) 
is more than offset by the increase in m (a factor of 7 r), 
and so the choice At = 7r is slightly more efficient (by 
roughly a factor of 2). Even larger At increase the range 
over which Cfe(At), rfc(At) vary, and hence require higher 
precision arithmetic (i.e. larger b). Consequently, it may 
not be practical to implement such large At. 

For the Lie-Trotter algorithm, we find that the error e 
is inversely proportional to to. As a specific comparison, 
to make e < 10 —5 , we needed p > 16 for the projection 
operator series, p > 11 for the reflection operator series 
(both with At = 7 r), and to > 2 21 for the Lie-Trotter 
algorithm. The computational cost 2 mpC of the series 
expansion algorithms is then of the order 7 x 10 2 C, which 
is a huge improvement over the corresponding cost mlC = 
4 x 10 6 C for the Lie-Trotter algorithm. The ratio of the 
two is consistent with the order of magnitude expectation 
(-eloge). 

To study the growth of the error with the evolution 
time, we varied the simulation time t , while holding At 
and p fixed. Our results are illustrated by FigEl For the 
series expansion algorithms, we find that e is proportional 
to t, implying that the errors of successive time steps ad- 
ditively accumulate, in accordance with Eg. (1311) . But we 
also find that for the Lie-Trotter algorithm such additive 
accumulation of error holds only for t < 1. Beyond that 
the error saturates with the saturation value proportional 
to At. This stoppage of error growth for large t indicates 
cancellations among the errors of different time steps, 
possibly due to correlations in the periodic evolution be¬ 
yond the first cycle (period of exp(— iHit) is 27 t) [37j. We 
note that for t > 1, we have roughly m = t/At oc t/e, 



FIG. 5: Dependence of the error e on the evolution time for 
the series expansion and the Lie-Trotter algorithms. The sym¬ 
bols □ , + and 0 respectively represent the results for the re¬ 
flection operator series (with p = 8, At = 1), the projection 
operator series (with p = 10, At = 1), and the Lie-Trotter 
algorithm (with At = 0.0001, 0.001 and .01 values connected 
by continuous, dashed and dotted lines respectively). 
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and not moc t 2 /e as per Eo. (l32l) . 


V. EFFICIENT SIMULATION OF LOCAL 
HAMILTONIAN EVOLUTION 

We now construct a rapidly converging series expan¬ 
sion for exp (—iHt). (The reason for decomposing the 
Hamiltonian into block-diagonal parts appears later in 
Section IV Cl ) It is well-known that an expansion in 
terms of the Chebyshev polynomials provides uniform 
approximation for any bounded function, with fast con¬ 
vergence of the series [38j. We use such an expansion 
for exp (—iHt), interpreting all matrix functions as their 
power series expansions [33l |. 


A. Chebyshev Expansion and its Complexity 

For any bounded Hamiltonian, its eigenvalue spectrum 
is within a range [A m i n ,A max ]. With a linear trans¬ 
formation, this range can be mapped to the interval 
[—1,1] that is the domain of the Chebyshev polynomi¬ 
als T n {x) = cos(ncos _1 x). Explicitly, 

g iHt _ g ^(•^max - l - -^min)^/2g iHt ^76) 

H = (2 H (A max A m i n )/)/(A max A m i n ),(77) 
t = f(A max - A min )/2 . (78) 

In situations where A m i n and A max are not exactly known, 
respectively lower and upper bounds for them can be 
used. Henceforth, we assume that such a mapping has 
been carried out and drop the tilde’s on H and t for 
simplicity. 

The Chebyshev expansion gives 

OO 

e ~ iHt =J2 Ck ^ W) ’ (79) 

k—0 

where the expansion coefficients are the Bessel functions: 

C 0 = - [ e~ ltcose dO = Jo(t) , (80) 

7T Jo 

C k = - [ e~ ltcose cos(ke) dO = 2(-i) k J k (t) .(81) 
7T Jo 

Note that the Chebyshev polynomials are bounded in 
[—1,1], and the coefficients Jk(t) = t k /(2 k k\) + ... fall off 
faster by a factor of 2 k compared to the corresponding 
coefficients t k /k\ of the Taylor series expansion. This is 
the well-known advantage of the Chebyshev expansion 
compared to other series expansions. 

For the special case analysed in Section IIVB II i.e. 
H = I — (Ri + R 2 )/ 2, Rf = I implies that the spec¬ 
trum of H is bounded in [—1,1], Then from the recursion 
relation for the Chebyshev polynomials, 


it follows that 

T k ((R1R2 -..)* + Wi • ■ ■)*) • 

(83) 

Thus we reproduce the series expansion obtained in 
Eo. (l63l) . with the same values for rk(t). Looking at it 
another way, the partial summation of the reflection oper¬ 
ator series for e~ lHt converts the Taylor series expansion 
into a better behaved Chebyshev expansion. 

When the Chebyshev expansion in Eq. (1791) is truncated 
at order p , its error analysis is identical to that in Section 
IIV B 21 Choosing t = mAt and At = 0(1), the constraint 
of Ea. dTUl) formally provides an efficiently converging se¬ 
ries with p = O(log(t/e)/log(log(t/e))). 

We have noted earlier that reflection operations are 
the largest evolution steps consistent with unitarity, and 
their use makes Grover’s algorithm optimal. Then, with 
e nzR /2 _ a g 00c i g ue ss for the evolution time step 
is At = 7 r. With this choice, the numerical tests in 
Section IIV Cl indicate that truncating the series at order 
p = 21n(t/e)/ln(ln(t/e)) is sufficiently accurate. 

A truncated series of the Chebyshev expansion is effi¬ 
ciently evaluated using Clenshaw’s algorithm, based on 
the recursion relation Eo. (15121) . One initialises the vec¬ 
tors \y P +i) = 0, | y p ) = C p \x), and then uses the reverse 
recursion 


I Vk) = C k \x) + 2 H \y k +i) - \yk +2 ) , (84) 

from k = p — 1 to k = 0. At the end, 


J2c k T k (H)\x) = (Col*) + \y 0 ) - |i/2»/2 (85) 

fc =1 


is obtained using p sparse matrix-vector products involv¬ 
ing H. The computational complexity of the total evo¬ 
lution is then 


0[mpCc) = O 


logft/e) 

log(log(i/e)) 



( 86 ) 


where Cc is the computational cost of implementing the 
recursion of Eo. (l84]) . 


B. An Alternate Strategy 


The Chebyshev expansion coefficients J k {t) are 
bounded for any value of f, unlike their Taylor series 
counterparts, and rapidly fall off for k > t. These prop¬ 
erties suggest an alternate evolution algorithm, i.e. eval¬ 
uate e~ lHt at one shot without subdividing the time in¬ 
terval into multiple steps [33 . 

The error due to truncating the Chebyshev expansion 
at order p is bounded by 


OO OO 

E i c *mi s Era 

k—p -\-1 k—p -\-1 


< 


t p+i 


( 87 ) 


t 


-1 


Tk+flH) = 2H T k (H) - T k _i(H) 


(82) 


2P{p + l)\ 


1 


2(p + 2) 
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provided the subleading contribution in Eg. (1681) can be 
ignored. The subleading contribution can certainly be 
neglected for p > f 2 /8, but the bound in Ea. ([571) may 
hold for even smaller values of p due to cancellations 
among subleading contributions of different terms in the 
series. Making Eq. (EZD smaller than e requires p > et/2, 
and the formal bound is p = 0(fe _1 A). The resultant 
computational complexity of the evolution, 

0 (pC c ) = Oite-^Cc) , ( 88 ) 

has a power-law dependence on e. But the power de¬ 
creases with increasing t, and so the computational com¬ 
plexity can be comparable to Eg. (1551) for values of e and 
t that are of practical interest. The extent to which 
the computational complexity would be enhanced by the 
need to control subleading contributions can be problem 
dependent, and needs to be determined numerically HTTl' . 

To implement this strategy, the Bessel functions Jk(t) 
upto order p need to be evaluated to 6 = H(log(p/e)) bit 
precision. That can be efficiently accomplished using the 
recursion relation, 

2 k 

= —Jk(t) — Jk+i(t) , (89) 

in descending order [isj. One starts with approximate 
guesses for J/(t) and J;+i(t), with l slightly larger than p, 
and uses the recursion relation repeatedly to reach Jq ( t). 
Then all the values are scaled to the correct normalisation 
by imposing the constraint Jo(t) + 2 ]>][■=? Jk{t) = 1- 
This procedure to determine the expansion coefficients 
requires O (pb 2 ) computational effort, and so does not 
alter the computational complexity. 

C. Digital State Implementation 

Summation of the series in Eq. (l79l) , truncated to order 
p, requires p executions of the Clenshaw recursion rela¬ 
tion, Eo. (1551) . Multiplication of a vector by a constant, 
and addition of two vectors, are easily carried out with 
the digital representation of Section fill El Multiplication 
of the sparse Hamiltonian with a vector, on the other 
hand, has to be carefully implemented such that quan¬ 
tum parallelism converts its computational complexity 
from classical O(N) to quantum 0(n). 

Multiplication by the diagonal elements of the Hamil¬ 
tonian has a trivial parallel structure, but its paral¬ 
lelisation for the off-diagonal elements of the Hamilto¬ 
nian needs decomposition of H into parts, with each 
part consisting of a large number of mutually indepen¬ 
dent blocks. As mentioned earlier in Section III A1 such 
a decomposition can be achieved for any sparse Hamil¬ 
tonian using an edge-colouring algorithm for the corre¬ 
sponding graph. With l colours, there are l Hamiltonian 
parts, each containing 0(N/2) mutually independent 2x2 
blocks. (Note that Hermiticity of the Hamiltonian re¬ 
lates the off-diagonal elements, Hjj +fl = HJ + •, that 


are represented by a single edge of the graph.) Evalu¬ 
ating the contribution of each Hamiltonian part in suc¬ 
cession, and combining the individual block calculations 
for each Hamiltonian part with a superposition of their 
block labels, the total computational effort for Hamil¬ 
tonian multiplication becomes 0(1 log(AT/2)) times the 
effort for a single 2x2 block multiplication. 

In the digital representation, the 2x2 block multiplica¬ 
tion becomes straightforward provided one can swap the 
6-bit register values, i.e. 

\j)\Vj) + 1 3 + M)l Vj+n) —* \j)\yj+u) + 1 3 + ■ (90) 

Were \yj) just a single bit, such a swap operation can be 
performed by the reflection operator, 

S = a x ®I , S 2 = I , (91) 

acting on the subspace {|j), \j + p)} ® {\yj), \Vj+n)}- For 
the 6-bit register, this swap operation has to be per¬ 
formed one by one for each bit. The swap can be easily 
undone after the off-diagonal element multiplication for 
a particular Hamiltonian part Hi, to use | y 3 ) again for 
the next Hamiltonian part. 

The digital circuit implementation of Eo. (l84l) . for a 
single Hamiltonian part Hi, is schematically illustrated 
in Fig® It has computational complexity 0(6 3 ) arising 
from evaluation of the Hamiltonian elements; the rest of 
the linear algebra operations have computational com¬ 
plexity 0(6 2 ). Including contributions of all the Hamil¬ 
tonian parts, and the computational effort needed to su¬ 
perpose the index j, we thus have the time complexity 

|o>6-^I (vk)j)b 

\c k ) b - 

I x j)b 

I*) -t-f- 

I j) -"- 

1 °) - 0 - 

I 0) b - \H} 

|0>6 - 

I ( yk+l)j)b 

\(yk+2)j)b - 

FIG. 6: Digital quantum logic circuit for executing the recur¬ 
sion relation of Clenshaw’s algorithm, Eq. (l84ll . to be executed 
with a uniform superposition over the index j. Operations for 
a single Hi containing only 2x2 blocks (labeled by j, j + /q) 

are shown. Among the controlled logic gates, H and de¬ 
note oracle operations specified by the Hamiltonian, is the 
swap operation of Eg. (1901) . fXl stands for the generalised Tof- 
foli gate implementing | a, b, c) —> |a, b,c+ ab ), and El labels 
the generalised C-not gate performing \a, b) —» |o, b — a). 


|(yfc+i)j+M>) 


- I j + im) 

- 1(2 Hi)^) 

K 2 

- \(Vk+i)j) 
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Cc = 0(lnb 3 ). We also point out that the space resources 
required to put together the full Chebyshev expansion are 
a fixed number of n-bit registers and O(p) 6-bit registers. 

Finally, note that the classical computational complex¬ 
ity for implementing Ea. (l84l) is Cc = 0(lNb 3 ). In our 
construction based on digital representation for the quan¬ 
tum states, the full quantum advantage that reduces N to 
n arises from a simple superposition of the quantum state 
label j, and this superposition in turn requires decompo¬ 
sition of the Hamiltonian into block-diagonal parts. 


VI. SUMMARY AND OUTLOOK 

We have presented efficient quantum Hamiltonian evo¬ 
lution algorithms belonging to the class P:P, for local effi¬ 
ciently computable Hamiltonians that can be mapped to 
graphs with bounded degree. Our construction exploits 
the fact that, the Lie-Trotter evolution formula can be 
reorganised in terms of reflection operators and Cheby¬ 
shev expansions (by partially summing up the BCH or 
the Taylor expansions), so as to be accurate for finite 
time step size At = 0(1). Specifically, P 2 = P and 
R 2 = I allow easy summation of a large number of terms, 
while the large spectral gap of P and R , due to only 
two distinct eigenvalues, provides a rapid convergence of 
the series [4(j ■ The net result is a dramatic exponential 
gain in the computational error complexity. Our expan¬ 
sions have better convergence properties than previous 
similar results [HJ, obtained by successively reducing 
the Hamiltonian evolution problem to simpler instances. 
Furthermore, our explicit constructions show how to de¬ 
sign practical efficient algorithms, and reveal the physical 
reasons underlying their efficiency. 

The formalism that we have developed has connec¬ 
tions to the familiar method for combining exponentials 
of operators, i.e. the Baker-Campbcll-Hausdorff formula. 
This formula can be partially summed up and simplified 
for exponentials of projection operators. Several identi¬ 
ties for projection operators that are useful in the process 
are described in the Appendix. In particular, the iden¬ 
tity of Eq. (1^1) may be useful in other applications of the 
BCH expansion. 

Our methods have introduced two concepts that go 
beyond the specific problem investigated here. One is 
that unitary time evolution using a large step size can be 
looked upon as simulation of an effective Hamiltonian. 
This effective Hamiltonian can be very different from the 
original Hamiltonian that defined the evolution problem 
in continuous time, as seen in our analysis of Grover’s 
algorithm. Such a correspondence between two distinct 
Hamiltonians that give the same finite time evolution is 
highly non-trivial, and underlies efficient summation of 
the BCH expansion. The technique of speeding up sim¬ 
ulations by finding appropriate equivalent Hamiltonians 
can be useful in a variety of problems defined as contin¬ 
uous time evolutions (including adiabatic ones). 

The other novel concept we have used is to map non¬ 


unitary linear algebra operations to unitary operators us¬ 
ing the digital state representation. High precision calcu¬ 
lations need a digital representation instead of an analog 
one. We have introduced such a representation for both 
the quantum states and operators, that maintains the ex¬ 
pectation values of all physical observables. It combines 
classical reversible logic with equally weighted linear su¬ 
perposition, and is essentially free of the unitarity con¬ 
straint for quantum states. Such digital implementations 
can help in construction of class P:P quantum algorithms 
for many linear algebra problems. 

A noteworthy feature of our algorithms is that they do 
not make direct use of any quantum property other than 
linear superposition—the constraint of unitary evolution 
is reduced to an overall normalisation that can be taken 
care of at the end of the computation and need not be ex¬ 
plicitly imposed at intermediate stages of the algorithm. 
Specifically, the digital representation of quantum states 
makes linear algebra operations involving action of block- 
diagonal Hamiltonians on a quantum state extremely 
simple. The Hamiltonian blocks can be processed in su¬ 
perposition on a quantum computer, while they can be 
handled by independent processors on a classical parallel 
computer. Consequently, our algorithms can be used for 
classical parallel computer simulations of quantum sys¬ 
tems, with the same exponential gain in temporal compu¬ 
tational error complexity. Of course, classical and quan¬ 
tum simulations will differ in the spatial resources, N 
classical variables vs. log(lV) quantum components, but 
the temporal cost will be identical with At = 0(1). In 
other words, classical and quantum complexities differ 
only in the cost C parametrising the resources required to 
carry out a sparse matrix-vector product, and the expo¬ 
nential gain in quantum spatial complexity simply arises 
when this product can be evaluated using superposition 
of an exponentially large number of blocks. 


Appendix: Some Identities for Projection Operators 

Let {Pi = \ei){ei\} be a set of normalised but not nec¬ 
essarily orthogonal projection operators: 

P 2 = Pl = P } , TriPiPj) = |(e i |ey>| 2 = |A^| 2 . (92) 

Functions of a single projection operator arc linear. For 
instance, the projection operators are easily exponenti¬ 
ated as 

exp (icfPi) = 1 + (e 1 ^ - 1 )P Z . (93) 

Furthermore, functions of two projection operators re¬ 
duce to quadratic forms (note that P t PjPi — |A ij\ 2 Pi). 
In general, the product of a string of projection operators 
reduces to an expression where each projection operator 
appears no more than once, because 

Pil Pi 2 Pil ■ ■ ■ Pi n Pil = ^ 11*2 ^* 2*3 • • ■ ■ ( 94 ) 
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Such simplifications reduce any series of projection oper¬ 
ators to finite polynomials, and various identities follow. 

(A) The operator (Pi — Pj) 2 has the orthogonal eigenvec¬ 
tors |ej)±|&j), with degenerate eigenvalue 1 — \Kj\ 2 - Also 
the operator [Pi,Pj] 2 has the same orthogonal eigenvec¬ 
tors, with degenerate eigenvalue |Ajj| 4 — |A^| 2 . These 
properties make both these operators proportional to 
identity in the subspace spanned by |e*) and | e j). In 
this subspace, therefore, we have (P,;. Pj} = Pi + Pj — 
(Pi — Pj) 2 = Pi + Pj — 1 -f- 1 Aiy| 2 , and the identity 


These expressions can also be derived using the identity 

[R,[R,(R,X]}) = [R,X] 

<+ (ad Pi) 3 (X) = ad P,(X) . (99) 

When A is also a projection operator, further simplifica¬ 
tion can be carried out using 

[Pi, [- Pi, P d ]] = Pi + Pj - 1 + |Ajj| 2 (l - 2 Pi), (100) 

in the subspace spanned by |ej) and \ej). 


e ±mPi e ±iirPj = (1 _ 2Pi)(l - 2Pj) 

= -l + 2\X ij \ 2 + 2[P i ,Pj\ 

= -exp ( ~ 2s hl (95) 

Eas. (|llll3l) correspond to the special case of this identity 
with A ij = 1 /VN. 

(B) In the subspace spanned by |e») and | ej), with a phase 
choice that makes A ij = (&i\e.j) real, | e*) ± | ej) are also 
the eigenvectors of Pi + Pj , with eigenvalues 1 ± A ij. The 
evolution taking | e*) to | ej) can therefore be achieved as 


(D) When A is a differential operator, interpreting 
df/dx = [d/dx,f], similar algebra yields 


p i<t>P, 


d(e I"_d_ -i4>Pi 

dx Idx' 


= -iPi^j- + (e - l)^p- + (2 - 2cos(j>)Pi^p- 
dx dx dx 


dcj) . . dPi 

= — iPi— - i sin 0—-1- (1 — cos 0) 

dx dx 


R 


dR 
dx . 


( 101 ) 


Note that P 2 = Pi leads to 


exp (-i(Pi + Pj)T)\ei) 

1 


2 L 
1 

2 


e -iT(l + Ay)(| e .) + |e . }) + e -iT(l-A y)( | e . ) _ | e . }) 

l e i) + l e j) + e 2iTXii (\ei) — \ej)) 

= -ie~ lT \ej), for T = 7r/(2Ay). (96) 


= i_g — iT(l+\ij) 


The Farhi-Gutmann search algorithm is the special case 
of this result with A ij = 1/VN- 

(C) The unitary transformation generated by a projec¬ 
tion operator is: 

e ,0g le -4P. 

= (1 + (e i4> - 1)P,)A(1 + (e"** ^ l)Pi) 

= X+i sin 0[P i; X] + (cos 0-1) ({P h X} - 2RXR) 
= X + i sin0[Pi, X] + (cos 0 — l)[Pj, [Pi, A]] . (97) 


In the Lie algebra language, the adjoint action of an oper¬ 
ator is defined as adY(A) = [Y, A]. So the above unitary 
transformation can be also expressed as 

e i<t>Pi Ae"^ Pi = e ^( ad P)(A) (98) 

= ^1 + i sin0(adPi) + (cos0— l)(adP;) 2 ^ (A). 


,, dPi dPi dPi 

p _ _ j_I p. — 

1 dx dx 1 dx 


R, 


R, 


dR 
dx . 


dR 

dx 


• ( 102 ) 


(E) The general BCH expansion for combining exponen¬ 
tials of non-commuting operators is an infinite series of 
nested commutators, and is cumbersome to write down 
at high orders. But it can be expressed in a compact form 
using exponentials of adjoint action of the operators (27j : 


e A e B = exp 


A + B 



(1-e 


ad AgS ad B\n 


n(n + 1) 



(103) 

When the operators involved are projection operators, 
the identities in (C), (D) convert exponentials of adjoint 
action of the operators to quadratic polynomials, effec¬ 
tively summing up infinite series. 

Similar simplification is possible for reflection oper¬ 
ators as well, due to the identity Ri = 1 — 2 Pj, e.g. 
(adPi) 3 (A) = 4 adPi(A). The resultant reorganised 
formula, if necessary with appropriate truncation, can 
be used as an efficient evolution operator replacing the 
Lie-Trotter formula. The series in Ea. (l53l) is an example 
of such a reorganised BCH expansion. 
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